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Monte Carlo radiation hydrodynamics: methods, tests and 
application to supernova Type la ejecta 



(N 

o 

(N 
W) 

< 

6 

(N 
> 

(N 
^sD 

^sD 
O 
(N 



13 



U. M. Noebauer, 1 * S. A. Sim, 2 M. Kromer, 1 F. K. Ropke 1 ' 3 and W. Hillebrandt 1 

1 Max-Planck-Institut fur Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany 

2 Research School of Astronomy & Astrophysics, Mount Stromlo Observatory, Cotter Road, Weston ACT 2611, Australia 

3 Institutfiir Theoretische Physik und Astrophysik, Universitdt Wiirzburg, Emil-Fischer-Str. 31, D-97075 Wiirzburg, Germany 



! March 2013 



ABSTRACT 

In astrophysical systems, radiation-matter interactions are important in transferring energy 
and momentum between the radiation field and the surrounding material. This coupling often 
makes it necessary to consider the role of radiation when modelling the dynamics of astro- 
physical fluids. During the last few years, there have been rapid developments in the use 
of Monte Carlo methods for numerical radiative transfer simulations. Here, we present an 
approach to radiation hydrodynamics that is based on coupling Monte Carlo radiative trans- 
fer techniques with finite-volume hydrodynamical methods in an operator-split manner. In 
particular, we adopt an indivisible packet formalism to discretize the radiation field into an 
ensemble of Monte Carlo packets and employ volume-based estimators to reconstruct the ra- 
diation field characteristics. In this paper the numerical tools of this method are presented and 
their accuracy is verified in a series of test calculations. Finally, as a practical example, we 
use our approach to study the influence of the radiation-matter coupling on the homologous 
expansion phase and the bolometric light curve of Type la supernova explosions. 
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1 INTRODUCTION 

In studying astrophysical objects, a detailed understanding and de- 
scription of the radiation field is vital, particularly if synthetic ob- 
servables are to be computed for comparison with observations. 
Conceptually, the radiation field in a fluid is not independent of 
the fluid state and their co-evolution has to be described self- 
consistently within the framework of radiation hydrodynamics. De- 
pending on the dynamical importance of the radiation field and 
the strength of the radiation-matter coupling, different strategies 
can be followed. If the energy associated with the radiation field 
is negligible compared to the total energy content, a de-coupled 
approach may be followed. For example, such a method has been 
used for the determination of synthetic light curves and spectra for 
Type la supernova (SN la) explosions a r ound maximum light (e.g. 
iKasen et al]|2006l : lKromer & SimlfcoogklJack et alj|201 lh . In cases 
where the radiative terms are dynamically important, however, a 
fully de-coupled treatment of the radiation field is not possible. 
For such applications a variety of different techniques have been 
used to follow the co-evolution of the radiation-matter state. In 
optically thick environments, the radiation field is well-described 
by the diffusion approximation and its evolution in radiation- 
hydrodynamical simulat ions can be incorporated by usin g flux lim- 
ited diffusion methods dLevermore &Pomranin3ll98ll) . This nu- 
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merical prescription is used, for exa mple, in the modelling of radia- 
tion d ominated accretion discs (e.g. lTurner et al.l2 003; Hiros e et"al] 
2009). In the opposite case of low optical depth, the influence of 
the radiati on field may be t reated by including a radia tive cooling 
term (e.g. Townsend 2009; van Marie & Keppens 2011), as is of- 
ten done in studies of stellar winds (e.g.jCarcia-Segur a et alll996l : 
iMellema & Lundqvistll2002h . In the intermediate regimes between 
the two extremes, a full radiation-hydrodynamical description of 
the radiation-matter state is necessary, for example when account- 
ing for convective motions in studies of ste llar atmospheres (e.g. 
IStein & Nordlundlll998l: lAsplund et al Il2000h. shock breakouts in 



Blinnikov et alj|2000l : iHoflich & Schaefej [2009 ; 



supernovae (e 

iPiro et al l2010h or when studyin g interactions of stellar expl osions 
with circumstellar material (e.g. lFrver et al.C OlO: Kaser feOlOl) . 



In this paper we present the numerical methods and the appli- 
cation of a new approach to radiation hydrodynamics that is based 
on Monte Carlo radiative transfer techniques. A similar strategy ha s 
been pursued in the calculations presented in lKasen"et afll l201il) . 
Monte Carlo methods have already sho wn tremendous success in 
pure radiative transfer a pplications (e.g. [Fleck & Cumrningslll97ll: 
Abbott &Lucvl 1 19851: iMazzali & Lucvl 1 19931; lLong & Kniggej 



20021: ICarciofi & Biorkmanl l200d IKasen et alj l200d: iHarrie; 



201 lh . Within this probabilistic approach, complex radiation- 



matter interaction physics can be simulated and problems with ar- 
bitrary geometries can be addressed. Here we aim to extend the 
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Monte Carlo method to radiation-hydrodynamical calculations and 
explore its practicality for modern astrophysical applications. 

The focus of this paper is to present the theoretical and nu- 
merical foundations and to verify the operation of our Monte Carlo 
radiation-hydrodynamical method. We begin with a brief overview 
of the theoretical concepts that govern radiating flows in Section 
[2] which is followed by an extensive description of the numeri- 
cal methods of our approach in Section [3] The physical accuracy 
and the computational feasibility of the techniques presented here 
are assessed in Section [4] in which the results of a series of test 
calculations are described. As a first application of the method in 
astrophysical environments we report in Section [5]on our investi- 
gation of SNe la ejecta. In particular, we study the influence of the 
radiation-matter coupling on the ejecta structure and the resulting 
effects on the bolometric light curve during the homologous ex- 
pansion phase. We summarize our results and conclude in Section 
i 



2 THEORETICAL BACKGROUND 

To model environments in which a significant part of the total en- 
ergy is stored in the radiation field, one must deal with the cou- 
pled evolution of the matter state and the radiation field. The for- 
mer changes due to external forces, gradients of the thermodynamic 
variables and the radiation pressure acting on the matter. Such ra- 
diation pressure gradients are the consequence of anisotropies in 
the radiation field whose temporal evolution is driven by its in- 
teractions with the surrounding medium. Generally, these inter- 
actions are strongly dependent on the state of matter, i.e. on its 
density, temperature and composition. The theory of radiating flu- 
ids provides an adequate self-consistent description of the dynam- 
ical behaviour of the radiation-hydrodynamical state of the cou- 
pled radiation field-matter system. In the following section, we 
will give a brief outline of some important aspects of this the- 
ory. In-depth discussions can be found in standard textbooks, e.g. 
iMihalas & Mihalasl dl984l) . 

To describe the energy and the momentum of radiating fluids, 
the standard hydrodynamical equations expressing conservation of 
momentum and energy are extended by including additional source 
terms that account for the influence of the radiation field. In astro- 
physical environments, the physical viscosity is typically insignif- 
icant compared to the viscosity inherent to the numerical schemes 
used to model fluid flows. Consequently, the ideal Euler equations 
are commonly employed. Modified by the influence of the radiation 
field, these take the form 

P^u = f-VP + G, (1) 
p^- t e = u-f-V-(Pu)+cG°. (2) 
Mass conservation, expressed by the continuity equation 

+ p(V ■ u) = 0, (3) 

is not affected by the radiation field. Here, p, u, e and P denote the 
fluid density, velocity, total energy and thermodynamic pressure, 
respectively. Possible external forces are accounted for by the force 
density /. The radiation field acts as an additional energy and mo- 
mentum source in the form of the radiation 4-force components G° 
and G (see below). Note that we have formulated the equations in 



terms of substantial derivatives 

which capture changes in the co-moving fluid frame. 

The G-terms essentially describe momentum and energy 
flows caused by an imbalance of absorption and emission inter- 
actions between the radiation field and its surrounding medium. 
Quantitatively, these interactions are characterised by the opac- 
ity, x( x ! tj n i ")> an d the emissivity r/(x, t; n, v) of the medium. 
These material functions depend on the frequency (y) and the prop- 
agation direction (n) of the radiation and will in general vary with 
time (i) and position (x), since the radiation-matter interactions de- 
pend strongly on the fluid state. The radiation force components can 
be specified in terms of these material functions and the specific in- 
tensity, I(x, t; n, v): 

G° = - 



G 1 = 



Au J dfl [x(as, t; n, v)I(x, t; n, v) — T](x, t; n, u], 

(5) 

(6) 

These can be understood as the net absorbed or emitted energy and 
momentum, respectively. The temporal evolution of the radiation 
field itself is in turn driven by the interaction with the environment 



+ n ■ V J I(x, t; n, u) 



18_ 

cdt 



rj(x, t; n, v) — x{x, t; n, u)I(x, t; n, v). 



(7) 



The combination of Equations {T}, 10, I0, (0 and an equation of 
state, relating thermodynamic pressure with internal energy, pro- 
vides the full set of radiation-hydrodynamical equations that de- 
scribe a radiating fluid. In this formulation, Equations (|5J and ((SJl 
describe how the energy and momentum transfer is obtained from 
the temporal evolution of the radiation field. In the following, we 
present in detail the numerical approach we developed to solve the 
radiation-hydrodynamical problem formulated by this set of equa- 
tions. 



3 NUMERICAL METHODS 

To determine the state of a radiating fluid, Equations ([]}, (0, 
0, and the equation of state have to be solved simultane- 
ously. Key to our approac h is the applicat ion of a simple operator 
splitting scheme (see e.g. iLeVequel |2005 . for a detailed descrip- 
tion of the operator splitting technique). In this Godunov-splitting 
framework, the temporal evolution is determined progressively, ac- 
counting for the pure hydrodynamical effects and the influence of 
the radiation field independently and in sequence. Specifically, in 
each time step, a new radiation hydrodynamical state is found by 
first performing a fluid dynamical calculation neglecting all radia- 
tive influences. For this step we employ a finite-volume hydrody- 
namical scheme, namely th e piecewise parabolic method (PPM, 
IColella & Woodward! [l984). This is followed by a second step 
which only accounts for the influence of the terms in the equa- 
tions governing the evolution of the radiation field and the fluid- 
dynamics terms are neglected. For this part of the simulation, we 
carry out time-dependent Monte Carlo radiative transfer calcula- 
tions, which allow us to evaluate the radiation terms and use them 
to update the hydrodynamical state. 
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The following sub-sections describe our scheme and the in- 
volved computational methods. We begin with an outline of the 
hydrodynamics solver in Section 13,11 followed by a detailed pre- 
sentation of the Monte Carlo radiative transfer techniques in Sec- 
tions l3,2l to l3.1 II The final step of updating the fluid state to account 
for the influence of the radiation field is described in Section [3.12| 

3.1 Hydrodynamical Calculation 

The hydrodynamical sub-problem of the operator-split ap- 
proach is solved with th e piecewise parabolic method of 
IColella & Woodward! dl984l) . This higher-order Godunov scheme 
is based on reconstructing a continuous fluid state from a discrete 
representation on a computational grid by a series of parabolas. 
In the spirit of finite-volume approaches, Riemann problems are 
defined at the cell interfaces by discontinuities in the fluid prop- 
erties, which arise from integrating the reconstructed fluid state 
over the domains of dependence, i.e. the regions that can influ- 
ence the interfaces. The solutions to the Riemann problems deter- 
mine the flux through the interfaces. By balancing the resulting in- 
and outflows in the grid cells, the temporal evolution of the fluid 
state can be calculated for one time step. This second-order re- 
construction scheme provides high er accuracy over the traditional 
constant-reconstruction method of iGodunovl ( 119591) . The detailed 
imple mentation of PPM in our radiation hydrodynamics code fol- 
lows lEdelmannl feOlti) . After determining the new fluid state, a 
Monte Carlo simulation is started to address the evolution of the 
radiation field as the second half of the splitting scheme. 

3.2 Monte Carlo Techniques 

In our Monte Carlo approach the radiation field is discretized into a 
large number of Monte Carlo quanta, hereafter referred to as pack- 
ets. Each of them carries a fraction of the radiation field energy 
and is propagated through the medium. The propagation path is 
determined stochastically but in accordance with the transfer equa- 
tion Q. From the ensemble of packet trajectories, all relevant radi- 
ation field quantities can be reconstructed. 

With the increasing availability of computational resources, 
Monte Carlo techniques have become a very popular and reward- 
ing approach to radiative transfer problems. Among the numerous 
applications of Monte Carlo radiative transfer techniques in astro- 
phy sics are the calculation of mass-loss rates in hot star winds (see 
e.g. I Abbott & Lucvll 19851 ; iLucv & Abbottll 19931 : Ivink et alj|200cl ; 
Sundqv ist et al. 2010), light curves and spectra of SNe la (see e.g . 
Mazzali& Lucvll 19931 ; lLucy||2005l ; iKasen et alj|200d ; lsirn]|2007l ; 
Kromer & Simll2009l) . ioni zation structure and s ynthetic spectra of 



photoionized nebulae (e.g. [Ercolano et al. 2003), mass-outflows in 
cataclysmic variables (e.g. Long & Knigge 2002; Noebau er et al.l 
120101) or active galactic nuclei (e.g. ISimll2005h . Compared with 
classical ray-tracing techniques, the Monte Carlo approach has cer- 
tain advantages. Most important, from a physical viewpoint, is 
the ease with which complex scattering and absorption processes 
can be incorporated. Since all interactions with the surrounding 
medium can be directly simulated during the propagation of the 
Monte Carlo packets, even the most complex atomic pr ocesses can 
be included in the radiative transfer calculation (see e.g. lLucy|2002L 
I2003ll2005l) . These interactions are all simulated locally in the co- 
moving frame, making the Monte Carlo algorithm entirely inde- 
pendent of the large-scale properties of the simulation and read- 
ily applicable even to problems with arbitrarily complex multi- 
dimensional geometries. 



Apart from these physically motivated advantages, the Monte 
Carlo method also brings computational benefits. As the propa- 
gation of one packet is independent of the behaviour of all oth- 
ers, Monte Carlo radiative transfer calculations can be easily paral- 
lelised and scale very well to large numbers of computational cores. 
This is of great significance since the efficient use of high perfor- 
mance computing facilities is an important consideration for the 
feasibility of modelling complex astrophysical systems. 

Of course, Monte Carlo radiative transfer methods also have 
their drawbacks. The accuracy and computational efficiency of the 
Monte Carlo approach are limited by the number of packets that 
discretize the radiation field and by the number of physical interac- 
tions they simulate. Consequently, whether Monte Carlo radiative 
transfer methods are appropriate s trongly depends on the specific 
problem under consideration (e.g. IPincus & Evans! 120091) . For ex- 
ample, in some applications, a detailed radiative transfer treatment 
is not required and the dynamical behaviour of the radiation field 
can be adequately addressed with approximate methods which per - 
form faster than the Monte Carlo approach (e.g. lKuiper et al 
Independently of the specific application, Monte Carlo methods al- 
ways introduce a certain level of statistical fluctuations in the sim- 
ulations. We shall return to the subject of minimising the influence 
of this Monte Carlo noise in Sections [3.10l andl 3.1 II 

3.3 Discretization 

As mentioned above, in the Monte Carlo approach the radiation 
field is discretized into packets. In early uses of the Monte Carlo 
machinery, e.g. in lAverv & House! dl968l) . the number of photons 
described by a packet was held const a nt thro ughout the simulation. 
However, we follow Abbo tt & Lucvl d 19851) and instead choose to 
discretize into packets of constant radiative energy. At every instant 
in the simulation, each packet represents a monochromatic parcel 
of radiative energy - i.e. a number of identical photons, all with 
a certain frequency v. When packets interact with the fluid, both 
the number of photons and the photon frequency associated with a 
packet can change but the energy it carries in the local rest frame 
of the fluid remains fixed. In addition, the energy packets are in- 
divisible - i.e. processes may create or destroy packets but never 
cause them to be split into multiple packets (see lLucvll2002L I2003L 
2005). This approach is motivated by the fact that total energy is 
conserved in interactions between the radiation field and its sur- 
rounding medium but the number of photons generally is not. The 
indivisible energy-packet method has been shown to be extremely 
powerful in solving radiative equilibrium problems owing to the 
ease with which it can ensure energy conservation (see e.g. ILucv! 
1999). By extending the method to include a net energy exchange 
between the radiation field and the matter, we continue to exploit 
the efficiency of this approach together with the properties of PPM 
to ensure global energy conservation in our simulations. The use 
of indivisible packets also avoids computational difficulties arising 
in cases where a fixed photon-number approach would require that 
packets are split (e.g. in modelling fluorescence or recombination 
cascades, where a physical process excited by a single photon leads 
to the re-emission of many). With the indivisible packet method, 
these processes are simulated with a probabilistic approach such 
that no packet splitting is required but that all the cascade channels 
are correctly sampled when a sufficiently large number of Monte 
Carlo packets are included (see Lucv 2 0021) . 

With this discretization scheme, the Monte Carlo packets are 
naturally well-suited to represent the mean intensity of the radia- 
tion field and its temporal evolution. In addition, all derived radia- 
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tion field characteristics can be easily formulated and reconstructed 
from Monte Carlo estimators and fully frequency-dependent opaci- 
ties could be readily implemented. However, the radiative flux is in 
general less accurately captured by the ensemble of energy packets. 
We will discuss the implication of this in more detail when consid- 
ering statistical noise and the accurarcy of our approach in Section 

urn 



3.4 Reference Frames 

An accurate and detailed description of the dynamical behaviour 
of the radiation field requires that special relativistic effects are 
considered. We have therefore designed our Monte Carlo radiative 
transfer algorithm to account for at least all first order special rela- 
tivistic effects. In principle, there are no obstacles to extending the 
method to higher order corrections. 

For handling relativistic effects, we must clearly distinguish 
between two reference frames. We define the spatial and tempo- 
ral discretization of the problem (i.e. our numerical grid) in the 
"lab" (or "observer") frame. The initialisation and propagation of 
the Monte Carlo packets is also performed in that frame. However, 
the natural choice for the treatment of matter-radiation interactions 
is the local rest frame of the fluid, which we will refer to as the "co- 
moving" frame. Henceforth, a subscript will be used to identify 
quantities that are defined in the local co-moving frame. 



3.5 Simplifications Adopted 

For the sake of clarity, in this section and below, we adopt some 
simplifications that apply to our current implementation and the 
test problems that we address in Sections 4 and 5. In particular, we 
restrict ourselves to one-dimensional problems, e.g. plane-parallel 
or spherically-symmetric media. By arranging the problem setup 
such that the symmetry axis coincides with the coordinate z-axis 
the radiation propagation direction can be specified by the scalar 



\x = n ■ e 2 



(8) 



In spherically symmetric geometries, fi is measured with respect to 
the radial direction. As a consequence of the symmetry properties 
of the problems we consider, all radiation-hydrodynamical quanti- 
ties only vary along one spatial coordinate. Thus, the components 
of the radiation force in the other two orthogonal directions (G 1 
and G 2 ) vanish and will not be considered further. 

In addition to the geometry restrictions, we treat radiative 
transfer in a grey, i.e. frequency-independent, approximation. Scat- 
tering interactions between the radiation field and the surrounding 
material are assumed to be coherent and isotropic. We stress, how- 
ever, that all these simplifications are not necessary - the approach 
can be readily generalized to multiple dimensions and frequency- 
dependent opacities. Finally, we use an ideal gas equation of state 
to relate fluid internal energy and thermodynamic pressure. 



3.6 Packet Initialisation 

At the beginning of a simulation, we need to generate an initial 
population of Monte Carlo packets that describe the initial radiation 
field. This generation process includes assigning each packet an 
initial position, direction and frequency. All the steps involved in 
this process have to accommodate the probabilistic nature of the 
Monte Carlo machinery. 

To initialize the population of Monte Carlo packets, an initial 



condition for the radiation field has to be chosen. As an illustrative 
example, we assume that the simulation is to be initialized in lo- 
cal thermodynamic equilibrium (LTE). In this case, the radiation 
energy density in the co-moving frame, Eo, follows the Stefan- 
Boltzmann law 



Eo = a R T 4 



(9) 



To correctly initialize the Monte Carlo packets, the total energy is 
first transformed into the lab-frame and summed over the entire 
computational domain 



Stot — ^ ' 



7 



1 + i/3 2 ) a R T 4 AV .. 



(10) 



Here, /3 and 7 denote the usual parameters of special relativity 

= -, (ID 

c 

1 



7 : 



(12) 



and the index i runs over all grid cells. The total energy £ tot is then 
divided equally into a chosen number of Monte Carlo packets (all 
packets are assigned the same initial lab-frame energy), which are 
spread over the grid cells according to the local radiative energy 
content. The initial position of a packet within a grid cell is chosen 
randomly. In LTE, the radiation field is isotropic in the co-moving 
frame, i.e. it has no angular dependence. Due to angle aberration 
effects, this isotropy is lost during the transformation into the lab 
frame. Consequently, the assignment of the initial propagation di- 
rection has to account for the angular dependence of the radiation 
field in the lab frame. In the grey approximation and under the re- 
striction to one-dimensional problems, the LTE lab frame specific 
intensity follows 



B(T) 



(13) 



7 4 (l-/3/*) 4 ' 

with B(T) denoting the frequency-integrated Planck function 

B(T) = — T 4 . (14) 

This angular dependence can be translated into a probability den- 
sity 

o2\3 -1 



■P 1 



2(l + l/3/3 2 ) (1-/3m) 4 ' 



(15) 



which can be sampled to give the relativistically correct directional 
distribution of the initial Monte Carlo packets. In the classical limit 
(/? — > 0, 7 — > 1), the density simplifies to 



(16) 



which can be easily sampled by the random number experiment 



/* = 2£ - 1, 



(17) 



where £ denotes a random number drawn uniformly from the inter- 
val [0, 1]. In the more general case, Equation H5\ must be sampled, 
leading to a more complex expression for /1 in terms of £, but the 
same principle applies. 

In general, each packet also has to be assigned a photon fre- 
quency ensuring that the packets represent the correct spectrum of 
the radiation field. In this work, however, this step can be skipped 
since we currently adopt a grey-approximation. A possible r ealisa- 
tion o f the more general sampling process can be found in iLucvl 
dl999l) . 
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Figure 1. Flow chart of the operator split algorithm and a detailed outline 
of the Monte Carlo radiative transfer step. During this process, the trajecto- 
ries of all packets are determined and the radiation field characteristics are 
calculated by the Monte Carlo estimators. 



3.7 Sequence of Monte Carlo simulations 

After their initialization at the beginning of the simulation, the 
Monte Carlo packets are propagated through the medium. During 
each time step, the packets are able to move (see Section [3~!8t . some 
packets are destroyed and others are created (see Section [3T9}. At 
the end of each time step, the properties of the currently active 
packets are stored so that these packets can be reactivated at the 
start of the next time step, after the fluid properties have been up- 
dated (see Section [3.12b . A graphical outline of the program flow is 
shown in Fig.Q] 



3.8 Packet Propagation 

In the Monte Carlo simulation, the packets propagate with the 
speed of light c through the medium. To simulate the dynamical 
evolution of the radiation field, the packets undergo interactions 
with the surrounding material as they propagate. In the Monte Carlo 
method, the propagation and interactions are treated stochastically. 
In particular, the location of packet interactions are determined 
by random number experiments - the trajectory of a propagating 
packet is terminated by an interaction with the surrounding medium 



when it covered a path-length 

X 

which follows from sampling the extinction law 



7 e 



-ix 



(18) 



(19) 



In addition to simulating physical interactions in this stochas- 
tic manner, two additional events that stem from the spatial and 
temporal discretization have to be taken into account. During the 
propagation process, packets can cross grid cell boundaries or reach 
the end of the current simulation time step. In the case of cell cross- 
ings, the changing fluid properties have to be taken into account, 
i.e. the interaction length I has to be recalculated. Since the Monte 
Carlo packets propagate with the speed of light, they will at max- 
imum travel a distance of d = cAt during a simulation cycle of 
duration At. At this point the propagation of a packet is suspended 
and its properties are stored in order to resume the propagation at 
the beginning of the next time step. 

3.9 Interaction Formalism 

Since our Monte Carlo packets represent photon packets, their in- 
teractions should model the physical interactions of photons with 
the surrounding medium: scattering, absorption and emission pro- 
cesses. When a packet experiences a physical interaction (i.e. once 
it has propagated the path-length given by Equation 1 18t. we must 
first determine the nature of the interaction. If we include both scat- 
tering (a) and absoiption («) contributions to the opacity 



X = o + k, 



(20) 



the relative strength of these will determine the probability of the 
interaction having been a scattering or absorption event. In particu- 
lar, we identify that a packet scatters if 



a + k 



(21) 



is fulfilled. To treat a coherent scattering event, the packet is trans- 
formed into the co-moving frame 



eo = £7(1 - /3/Li) 



(22) 



where the packet energy is conserved and a new propagation di- 
rection fi' is drawn isotropically (cf. Equationll7). Afterwards, the 
packet properties are transformed back into the lab frame 



e = £o7(l- + /3/Uo) 



(23) 
(24) 



1 + /Vo ' 

to resume the propagation. For all other outcomes of the experiment 
(Equationl21t. an absorption event occurs, resulting in the destruc- 
tion of the packet. Consequently, this packet stops its propagation 
and is no longer considered in the remaining simulation process. 

Destruction of packets by absorption interactions drains the 
energy content of the radiation field. However, radiation energy is 
also created due to emission by the medium. In the case of thermal 
emission, the radiative energy pool in a grid cell of volume AV is 
increased by 



AS = 47K<T R T 4 AVAt 



(25) 



during a time step of size At assuming constant temperature and 
opacity. To represent this effect in the Monte Carlo calculation, 
a number of new packets are created and launched during each 
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cycle whose total energy is consistent with this energy injection. 
The initial packet properties (e.g. direction of propagation, which 
we assume to be isotropic in the co-moving frame) are assigned 
by sampling appropriate probability distributions, in analogy to the 
process of representing the initial radiation field at the onset of the 
simulation (see Section [3T6l l. However, these emitted packets are 
not all injected into the simulation at the start of a time step but at 
randomly determined times during the time step, thereby account- 
ing for the continuous energy injection by the thermal emission 
process. 

In all interaction processes, momentum and/or energy are 
transferred between the radiation field and the surrounding mate- 
rial. Through their interaction behaviour, the packets account di- 
rectly for the impact of the momentum and energy flows on the 
radiation field. The complementing effect on the surrounding fluid 
material is equally important (see Equations [I] and [2). These flow 
terms can be reconstructed from the Monte Carlo simulation with 
the help of volume-based Monte Carlo estimators (see next Sec- 
tion, [013. 



3.10 Monte Carlo Estimators 

The Monte Carlo packets were introduced to discretize the radi- 
ation field and, during their propagation, simulate its interactions 
with the medium. In principle, all properties of the radiation field 
can be directly determined from the instantaneous properties of the 
packets at any given time, e.g. from their instantaneous distribu- 
tion in real space, momentum space and frequency space at the end 
of a time step. However, the accuracy with which radiation quan- 
tities can be reconstructed is limited by Monte Carlo noise and it 
is advantageous to consider how t o extract max imum information 
from the Monte Carlo simulation. iLucvl i 19991) showed that radi- 
ation field quantities can be more accurately determined by using 
volume-based Monte Carlo estimators rather than directly counting 
properties of the packets. In this approach, the complete ensem- 
ble of trajectories that packets traverse during a time step is used 
to reconstru ct the radiation field characteristics (see lLucy|[l999l . 
120031 120051) . Thus, each packet contributes to the cell-averaged 
values of the radiation field characteristics according to the frac- 
tion of the time that it spent in the cell. Compared to reconstruction 
methods based on considering the final packet distribution at the 
end of the propagation step, this cumulative approach significantly 
reduces the Monte Carlo noise, which is in general the limiting 
factor in the applicability of Monte Carlo methods. To reconstruct 
the cell-averaged radiation field energy density all packet trajec- 
tory segments that cross the cell under consideration or lie within 
it are taken into account. On each of these trajectory segments, 
the packet contribution scales with its energy e. Every contribu- 
tion is weighted according to the ratio of the time 5t% the packet 
spent on the trajectory element to the full duration of the simula- 
tion time step At. Replacing the propagation time with the path 
length U. = cSti of the individual segments and summing up, gives 
the estimator 



E = 



AVcAt 



5> 



(26) 



for the radiation field energy density. Here, the summation runs 
over all trajectory elements in one cell, implying that a packet may 
contribute several times to the estimator. Using the fundamental 
relation between the radiation field energy density and the mean 



intensity 
J 



(27) 



a Mo nte Carlo estimator for the latter can be formulated (see lLucvl 
I1999L equation 12) 



3 : 



1 



■5> 



(28) 



AitAVAt 

By restricting the estimator sum to contributions by packets propa 
gating into a certain directional bin [fi, n + A(m], 

1 



27rAVAt 



[m>M+ a /j] 



el. 



(29) 



the specific inte nsity can also be determined in a similar manner 
(see lLucvl2005L equation 2). With this expression, estimators for all 
radiation field characteristics that depend on the specific intensity 
or its moments can be easily formulated. 

For radiation-hydrodynamical problems, the radiation force 
components have to be reconstructed from the Monte Carlo sim- 
ulation to determine the energy and momentum flow between the 
fluid and radiation field. The components of the radiation force can 
be interpreted as the difference between absorbed and emitted radi- 
ation energy and momentum respectively. We clearly separate the 
contributions to the radiation force terms caused by scattering in- 
teractions from the energy and the momentum transferred in pure 
absorption and emission events. To identify the latter contribution, 
each packet trajectory is considered to affect the cell-averaged ab- 
sorbed energy, even if the packet did not explicitly interact dur- 
ing the propagation cycle. Conceptually, we are therefore count- 
ing all absorption events that could have happened in the simula- 
tion (weighted by their probability of occurring), rather than sim- 
ply counting the events that did happen. Thus, while propagating a 
path-length a packet contributes with 

(30) 



and with 

[line 



Ap 



cAV 



(3D 



to the cell-averaged absorbed energy and momentum densities re- 
spectively. The complementing energy and momentum injection 
due to thermal emission is determined analytically for each cell. 
Accounting for the frame transformation into the lab frame the in- 
jection rates can be formulated as 

E = 47k<t r T 4 , (32) 

p = i/3 7 Kcr R r 4 . (33) 

c 

In addition to determining the contribution of pure absorption 
and re-emission events to the energy and momentum transfer terms, 
scattering interactions have to be incorporated. We consider scatter- 
ing interactions formally as a combination of an absorption event, 
followed immediately by the coherent and isotropic re-emission of 
a Monte Carlo packet. This interpretation allows us to reconstruct 
the scattering contribution analogously to Equations ( 130) and J3 1 b - 
In general, however, formulating an analytic expression to quan- 
tify the re-emission part is non-trivial. This difficulty can be cir- 
cumvented by exploiting the fact that Monte Carlo packets are not 
destroyed in scattering interactions. Thus, the packets can simulate 
both the absorption and re-emission parts of the scattering events 
at the same time. By additionally assigning each packet a set of 
hypothetical "post-scatter" properties, i.e. drawing a propagation 
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direction [i' that translates into an energy e' and a direction fjf in 
the lab frame according to the transformation rules (1231 and ( 124b . 
each packet trajectory I contributes to cell averaged energy and mo- 
mentum transfer according to 

ol 



Ap = 



cAV 



{[IE 1 - [l 1 E 1 ). 



(34) 
(35) 



Here, the superscripts i and / denote the current packet properties 
during the propagation and the formal post-scattering state, respec- 
tively. 

Finally, volume-based Monte Carlo estimators for the radia- 
tive force terms G° and G 3 are obtained by gathering all contribu- 
tions of scattering and pure absorption/emission events 



G° = 



AVcAt 



G 



AVcAt 



(36) 



(37) 



Once transformed into the local co-moving frame, these estimators 
take the expected form for the corresponding radiation field char- 
acteristics measured in this frame. In particular, the estimators ((36) 
and l !37t become 



Go 



> KqIoEo KoB(T), 

L ^ p 



r 3 - 

(_r — 



AVocAto 
1 



(38) 
(39) 



AVoeAfo ^— ' 

after transformation into the local co-moving framed using 

Go = 7 (G°-/3G 3 ), (40) 
Gl = 7 (G 3 -/?G°). (41) 

These co-moving frame estimators for G[] and Go reproduce ex- 
actly the radiative energy and momentum sources in this frame, 
which can be formulated as 

4tt 



G°o = —K (Jo-B(T)) 
G 3 = *°F . 



(42) 
(43) 



To obtain the same transformation behaviour for the Monte 
Carlo estimators in spherically-symmetric geometries, the direction 
scalar has to be averaged over the trajectory segment. 

3.11 Monte Carlo Noise 

Due to the stochastic nature of the Monte Carlo technique, all as- 
sociated quantities are subject to statistical fluctuations, limiting 
the accuracy of the radiative transfer calculation. The obvious ap- 
proach to reduce the Monte Carlo noise by increasing the number 
of active packets is limited by computational efficiency consider- 
ations. In this context, the concept of volume-based Monte Carlo 
estimators, as described in the previous section is vital since it al- 
lows us to extract the maximum amount of information from the 
propagation behaviour of a given number of radiation packets. The 
reduction of the statistical fluctuations in the reconstructed quanti- 
ties is significant but still limited. In particular, the noise level in 
the energy-momentum transfer and thus the accuracy with which 



For detailed derivations of t he transformation la ws for the radiation field 
characteristics, see e.g. Mihalas & Mihalal fl984l) . 



the radiation-hydrodynamical coupling is captured, varies with the 
state of the radiation field, due to our discretization scheme and 
to the form of the radiation force components in the co-moving 
frame. Reconstructing the energy transfer from a difference of two 
terms as in Equation d42t . yields accurate results if the contribu- 
tions are clearly separable. However, if the mean intensity deviates 
from the equilibrium configuration (i.e. B(T)) only on the level of 
the Monte Carlo noise, the resulting heating/cooling term will be 
obscured by statistical fluctuations. Similarly, the momentum trans- 
fer is only accurately captured if the radiative flux is well-resolved 
(cf. Equation I43t, To achieve high accuracy in the reconstruction 
via the estimator formalism, the packet contributions to the estima- 
tors have to be much smaller than the reconstructed quantity. This is 
the case for the mean intensity, since our discretization scheme (in- 
divisible energy packets) is most naturally suited to well-represent 
this radiation field characteristic. In the reconstruction of the ra- 
diative flux, however, the packet contributions to the estimator are 
weighted with the propagation direction. As a consequence, the ra- 
diative flux is only well-resolved in our approach if Fo < 4-7T Jo, 
but obscured by Monte Carlo noise for Fo <C 47r Jo, i.e. for nearly 
isotropic radiation field configurations. 

In summary, our Monte Carlo approach is ideal for describ- 
ing the mean intensity and our estimators will accurately capture 
the radiation forces when significant deviations from LTE or from 
isotropy are present. However, in cases where the radiation field is 
nearly isotropic or close to the LTE configuration, the radiation hy- 
drodynamical effects, as described in the radiation force terms, may 
be subject to significant statistical fluctuations. In order to achieve 
a further suppression of the noise in this regime, we have found it 
useful to smooth the reconstructed radiation field quantities over 
neighbouring grid cells. This additional noise reduction method 
will be of importance when addressing SN la ejecta (see Section 
[53}. 



3.12 Dynamical Influence of the Radiation Field 

With the numerical methods described in the previous sections, the 
temporal evolution of the radiation field is determined in a Monte 
Carlo step as part of the operator splitting approach to radiation 
hydrodynamics. After solving the fluid dynamical evolution using 
PPM and simulating the radiative transfer in the Monte Carlo sim- 
ulation, the calculation of the new radiation-hydrodynamical state 
ends with the inclusion of the radiative influence on the dynam- 
ics specified by G° and G 3 . These terms are reconstructed from 
the Monte Carlo simulation using the estimators defined in Equa- 
tions d36b and i37\ and describe the energy and momentum transfer 
onto the surrounding material. Consequently the fluid energy and 
momentum are updated according to 



p m u=G ' 

P Wt e = cG , 



(44) 
(45) 



as the second and last step of the Godunov splitting approach. This 
concludes the calculation of one time step and a new simulation 
cycle is entered by stepping through the segments of the splitting 
scheme again, starting by solving the fluid dynamics for the next 
time step (see Fig.[T}. 
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4 TESTING 



The methods described in the previous section have been imple- 
mented into a numerical code, hereafter referred to as MCRH. Cur- 
rently, the program is able to calculate the temporal evolution of 
radiating fluids in an Eulerian or Lagrangian reference frame, as- 
suming grey radiative transfer and one-dimensional geometries (see 
Section [3~5t , Before using the code to model radiative flows in as- 
trophysical environments or implementing more complex geome- 
tries and opacities, the computational feasibility and physical ac- 
curacy of our approach and its implementation has to be veri- 
fied. PPM is commonly used to simulate fluid flows in astrophys- 
ical applications, e.g. in Type la and II supernova explosions (e.g. 
Ropk e & Hillebranddl2005l:|janka & Muellejl 19961) or in relativis- 
tic jets (e.g. iMarti et alJI 19971) and has already been tested exten- 
sively. Consequently, the main focus of this discussion lies in test- 
ing the radiative transfer and the radiation-matter coupling compo- 
nents of our approach. 

We start this verification process by first testing the propaga- 
tion behaviour of the Monte Carlo packets. In the next stage, the 
interaction mechanism is checked by examining the equilibration 
behaviour in a series of toy simulations. Finally, we turn towards 
radiative shocks, the standard test problem for radiation hydrody- 
namics. In calculating sub- and supercritical shocks, the workings 
of both the radiation-matter coupling and our complete radiation 
hydrodynamics method can be verified. 

4.1 Random Walk 
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Figure 2. Comparison between the theoretically expected diffusion be- 
haviour (solid lines) of a centrally peaked radiation field in an optically 
thick, pure scattering medium and the corresponding Monte Carlo sim- 
ulation (circles). The lower panel shows the absolute difference between 
the theoretocally expected (-E"th) and simulated (E s \ m ) energy densities, 
A-Er = i? t h — -B s i m . 



results are compared with the theoretically expected behaviour set 
by Equation l |47t . The excellent agreement between the simulation 
and the theoretical predictions demonstrates the accurate operation 
of the basic Monte Carlo processes driving the packet propagation. 



As presented in Section [3751 the Monte Carlo packets undergo a 
multitude of physical interactions and numerical events on their 
propagation trajectories. Since this behaviour is essential for accu- 
rately simulating the temporal evolution of the radiation field, our 
first test calculation aims to verify the pack et propagation process. 
For this purpose, we follow iHarriesI (120 111) and consider a purely 
scattering medium in the optically thick limit. The entire radiation 
field energy is initially concentrated in the central cell of the com- 
putational domain (which has length L), mimicking a 5-function. 
In this problem, the initial radiation peak disperses according to the 
diffusion equation 



AE_ _ 

At 

yielding 



D 



d 2 E 
dz 2 ' 



E(z,t) = E(t ) 



V^nDt 



exp 



--) 
4DtJ 



(46) 



(47) 



With interactions restricted to scattering only, the diffusion coeffi- 
cient is given by 



(48) 



In a Monte Carlo simulation of this problem, the packets are ex- 
pected to perform a random walk due to the large scattering op- 
tical depth. Consequently, the temporal evolution of the radiative 
energy constitutes a suitable test to verify the correct propagation 
behaviour of the Monte Carlo packets. 

Following lHarriesI d201 ll) . the test simulation is performed in a 
planar symmetric box of length L = 1 cm that is divided into 101 
equally spaced cells. In the central cell, an initial radiation field 
energy density of E = 10 10 erg cm -3 is deposited and discretized 
by 10 5 Monte Carlo packets. Fig. [2] shows the temporal evolution 
of the radiative energy density in our Monte Carlo simulation. The 



4.2 Equilibration Behaviour 

Equally important as the correct propagation behaviour is the accu- 
rate calculation of the momentum and energy transfer between the 
radiation field and the surrounding material. We verify this part of 
our scheme by examining the relaxation behavio ur towards equi- 
librium in a series of toy simulations. Following iTurner & Stonel 
j200lh and lHarriesl ( l20T l|), we consider a radiating fluid initially far 
from equilibrium. In these particular tests, only absorption and ther- 
mal emission events occur in the medium and all influences of the 
thermodynamic and radiation pressure can be neglected. Therefore, 
the dynamical behaviour of the radiating fluid can be expressed in 
terms of the radiation field energy density _Er an d the internal en- 
ergy density Eg, which follow (see Harries 201 1, equation 23) 



^-E G = ck£ r - Akcj^T^Eg) 
at 

and 



(49) 



(50) 



As a first test, we follow the relaxation behaviour of the 
radiation field and the internal gas energy in a one-dimensional 
plane-parallel box of length L = 4 cm, which is discretized on 



4 cello F° r this, we adopt the parameters p = 10 7 gem 3 , 
E R (0) = 10 12 erg cm -3 , k = 4x 10~ 8 cm -1 ,u = Oandtheadi- 
abatic index 7 = 5/3. For these parameters, the equilibrium value 
for the internal energy of the fluid is E G q = 4.2 x 10 7 erg cm -3 . 



2 The finite-volume hydrodynamical solver based on PPM requires a do- 
main size of at least 4 cells due to the pa rabolic reconstruction and the 
discontinuity detection ( Colella & Wood ward 19841) . 
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Figure 3. Results of the equilibration simulations described in Section I4T21 
The blue symbols show the results of the Monte Carlo simulation with an 
initial internal energy density of Eq(0) = 10 2 erg cm -3 for a series of 
selected snapshots. The corresponding results obtained for the initial energy 
density of Eq(0) = 10 10 erg cm -3 are shown with red circles. The tem- 
poral evolution of the internal energy density as predicted by Equation (49) 
for both cases is indicated by the solid black lines. The initial radiation field 
energy density is set to i?n(0) = 10 12 erg cm -3 in both simulations. In 
the lower panel, the relative difference between the theoretically expected 
(i?th) and the simulated (i? s i m ) evolution of the internal energy density, 
(£ t h — E s i m ) I -Eth, is shown for both calculations. 



We perform the relaxation test twice, first with the initial internal 
energy of the fluid set to Eg(0) = 10 10 erg cm -3 and then with 
Eg(0) = 10 2 erg cm -3 . To predict the evolution of the systems, 
only Equation l |49t has to be considered, since the radiation field 
energy can be assumed to be constant in this configuration. The 
theoretical internal energy evolution, obtained from numerical in- 
tegration of Equation j49\ for both setups is compared with the 
results of the corresponding Monte Carlo simulations in Fig. [3] As 
illustrated in detail in the lower panel, the agreement is excellent. 
Only a systematic deviation persists, which is caused by the finite 
time resolution in the simulation. In the operator- split approach, 
the gas temperature is assumed to be constant during the radiation 
transfer cycle. Thus, the emissivity and in turn the amount of emit- 
ted energy is slightly overestimated. 

In a second test calculation, the initial radiation field en- 
ergy content was set to zero and the internal energy density to 
Eg(0) = 10 8 erg cm -3 . The remaining simulation parameters are 
adopted from the previous equilibration setup. Over time, the fluid 
will drive the radiation field towards equilibrium by the emission 
of thermal photons. To theoretically predict the evolution of the 
radiation-matter state, the coupled differential Equations i50\ and 
j49\ can be solved simultaneously by numerical integration. In con- 
trast to the previous simulations, neither the radiation temperature, 
nor the gas temperature can be assumed to be constant. In Fig. [4] 
the results of our Monte Carlo simulation are tested against the the- 
oretically predicted evolution. Again, only systematic deviations, 
which are caused by the finite duration of the simulation time steps, 
persist on the per cent level. 

The agreement between our test calculations and the theoreti- 
cally predicted behaviour for the relaxation towards equilibrium is 
excellent as illustrated in Figs.[3]and[4] indicating that our numer- 
ical methods describing the matter-radiation interactions and their 
implementation are accurate and work correctly. 



Figure 4. Equilibration calculation with an initially absent radiation field 
and a starting internal energy density of Eq = 10 8 ergcm~ 3 . For a se- 
ries of selected time snapshots, the theoretical evolution of this system, ob- 
tained by numerically integrating Equations i50] and i49\ is displayed as 
black lines. The corresponding results form our Monte Carlo simulation are 
depicted as blue and red circles for the radiative and internal energy density, 
respectively. As in Fig. [3] the relative difference between the theoretically 
expected and the simulated energy evolution is shown in the lower panel. 

4.3 Radiative Shocks 

The calculations presented above probed the propagation of the ra- 
diation field, the basic interaction machinery and the determina- 
tion of heating and cooling terms, neglecting any effects of radia- 
tion and thermodynamic pressure and the fluid motion. In our final, 
most challenging test, we verify that the interaction processes and 
our approach as a whole yield the correct radiation-hydrodynamical 
coupling. For this purpose we examine sub- and supercritical ra- 
diative shocks. The semi nal works describin g thi s type of shock s 
analyti cally date back to lZei'dovichl dl957ah and lRaizej Jl957alf 5 l 
and to Marshak (1958). The results of the two former studies are 
summarized in detail bv lZel'dovich & Raized dl967t) . 

Radiative shocks exhibit a structure and dynamical behaviour 
that is very different from their purely hydrodynamical counter- 
parts, due to the presence of a radiative precursor. The pre- shock 
material is heated by the flow of radiative energy through the shock 
front. Depending on the strength of the shock and in turn on the 
amount of heating of the pre-shock material, two classes of ra- 
diative shocks can be dis tinguished. Following the discussion of 
Miha las & Mihalaj dl984l) . we denote the gas temperature behind 
the shock front with T2 and immediately in front of it with T_ . In 
the case of T_ being much lower than T2, due to the small amount 
of pre-heating, a sub-critical shock is encountered. With increasing 
shock strength, equivalent to the shock front moving faster, the ra- 
diative precursor heats the material more and more until it reaches 
the critical configuration of T_ = T%. At this point, a further in- 
crease in the shock strength only results in a deeper penetration of 
the radiation precursor. Consequently, T_ never surpasses the tem- 
perature behind the front in such supercritical shocks. 

The structure of radiative shocks has already been de- 
termined in various nu meri cal studies, inclu ding works by 
iHeaslet & Baldwinl dl963h and lSincell etai] dl999h . Following the 

3 The corr esponding trans lations to English can be found in IZerdovict] 
ll957bl) and lRaizej jl957bl) . 
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Figure 5. Temperature profiles of a sub-critical radiative shock calculated 
with our Monte Carlo radiation hydrodynamics (red lines) program and the 
ZEUS-MP2 code (blue lines). The gas temperature is depicted with solid 
lines, while the radiation temperature is shown as dashed and dotted lines. 
Profiles are shown for t = 5.4 X 10 3 , 1.7 X 10 4 , 2.8 X 10 4 and 3.8 X 
10 4 s and are plotted in the rest frame of the un-shocked material. Note the 
Monte Carlo noise in the radiation temperature ahead of the shock. Here, 
the radiation field is only represented by a small number of packets. 

proposal of lEnsmanl (1994), these shocks have been used as a 
common test problem to verify the accuracy of numerical ap- 
proaches to radiation hydrodynamics. In particular, the realisa- 
tion of the matter-radiation coupling in th e ZEUS-code has been 
tested extensively using radiative shocks jTumer& Stone] 1200 ll ; 
lHaves & Norman 2003; Ha ves et al1l2006t) . For comparison to our 
Monte Carlo radiation-hydrodynamical simulations, we will use 
results calcu lated with the latest version of this code, ZEUS-MP2 
l lHaves et alj2006t) . 

For the calculation of sub- and supercritical radiative shocks 
we choose the p roperties of the radiating fluid according to 
lHaves et allfeOOd) which in tur n were motivated by the simulations 
performed by lEnsmanl d 19941) . In particular, a one-dimensional 
plane-parallel domain of length 7 x 10 10 cm and a fluid with ini- 
tial density of 7.78 x 10 _8 gcm -3 is considered. The tempera- 
ture is set to 85 K at the inner, reflecting boundary and decreased 
linearly to 10 K at the outer open boundary of the domain. We 
adopt this profile to f acilita te comparison with the calculations per- 
formed by Ensmanl l ll994l) . in which a temperature gradient had 
to be imposed for reasons of numerical stability. Throughout the 
domain, a uniform grey absorption opacity of 3.1 x 10 -10 cm -1 
is chosen, resulting in a photon mean free path that is roughly 20 
times shorter than the extent of the computational box. The shock 
is created by driving the material into the reflecting inner bound- 
ary with a bulk velocity of —6 km s~ 1 for the sub-critical case and 
with —20 kms -1 for the supercritical calculation, respectively. We 
determine the evolution of the resulting radiative shocks with our 
Monte Carlo program and as a reference with the ZEUS-MP2 code. 
Figs.[5]and|6]show comparisons between the gas and radiation tem- 
perature structures calculated with both codes for the sub - and the 
superc ritical shocks respectively. To be compatible with lEnsmanl 
d 19941) we display the shock structure in the rest fram e coordinates 
of the un-shocked material l lHaves & Normanl 120031 ; lHaves et al.l 
l2006h . 

As shown in Figs.[5]and[6] the radiative shocks calculated with 
our MCRH code exhibit the expected overall structure and charac- 




Figure 6. Analogous to Fig.|5]but now displaying the corresponding calcu- 
lation of the supercritical shock. The results from the Monte Carlo simula- 
tion are again shown in red and the ones obtained with ZEUS-MP2 in blue. 
Gas and radiation temperature are indicated with solid and dashed/dotted 
lines respectively and shown with respect to the un-shocked material. Tem- 
perature profiles are presented for t = 8.6 X 10 2 , 4.0 X 10 3 , 7.5 X 10 3 
and 1.3 X 10 4 s. 

teristic features. In the subcritical case, a weak radiative precursor 
leads to mild heating of the pre-shock material. The effect of the 
precursor is much more prominent in the second case of the su- 
percritical shock. Here, the sharp shock front is washed out by the 
strong heating of the upstream material by the radiative flux pene- 
trating the pre-shock domain. As predicted for supercritical shocks, 
the material immediately ahead of the shock is heated significantly, 
reaching the same temperature as behind the front. Overall, our nu- 
merical results agree very well with the ZEUS-MP2 calculations, 
especially with respect to the location of the shock and the gas 
temperature profiles. However, our Monte Carlo simulations pre- 
dict a stronger and deeper penetration of the unshocked material by 
the radiation field. This results in increased heating of the material 
and differences in the radiation temperature profiles with respect 
to the ZEUS-MP2 results. These differences are most likely caused 
by the simplifications implemented in the ZEUS-MP2 code, which 
treats the radiative flux in the diffusion a pproximation, augm ented 
by the introduction of a flux limiter (see lHaves et al.ll2006l equa- 
tion 1 1). The influence of this simplified sche me on the struct ure of 
radiative shocks has already been studied by lEnsmanl d 1994b . Cal- 
culating the evolution of a supercritical shock on the one hand by 
relying on the diffusion approximation and on the other hand by 
solving the radiation-hydrodynamical equations without any sim- 
plifications direc tly, revealed the same behaviour of the radiative 
precursor (see lEnsmanll 1 994j fig. 15). Despite its statistical nature, 
the Monte Carlo approach provides a direct solution to the trans- 
fer equation without introducing any physical simplifications. For 
these reasons we believe that the differences are due to the em- 
ployed methods and not a shortcoming of our approach to radiation 
hydrodynamics. 

The test calculations presented in this section have been car- 
ried out without using the smoothing capability of our program (see 
Section l3.1 U . In the simulations determining the structure of radia- 
tive shocks, the noise suppression of the volume-based estimator 
approach was sufficient to accurately resolve the heating effects in 
the precursor regime with large deviations from LTE. For example, 
in our simulation of the subcritical shock a maximum of 6 x 10 5 
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packets were active at a time, simulating the evolution of the radi- 
ation field. However, as anticipated in Sections [3.31 and [3 . 1 1 1 in the 
regions around the shock front, where the radiation field is close to 
LTE, the radiation force components are subject to non-negligible 
statistical fluctuations. But this noise component is effectively sup- 
pressed on the characteristic hydrodynamical time scales, due to a 
high packet recycling rate. Typically, the entire packet ensemble is 
re-populated multiple times during a radiative transfer step, causing 
a smooth temperature profile even in the near-LTE regions. 

With the calculations presented in this section, the differ- 
ent individual mechanisms of our method have been successfully 
tested and the correct radiation-hydrodynamical behaviour of the 
approach as a whole has been demonstrated without applying any 
smoothing. A first application to an astrophysical environment is 
presented in the next Section. 



5 APPLICATION 

Our framework has been specifically designed to address radiation- 
hydrodynamical problems in astrophysical systems with large out- 
flows such as supernova (SN) explosions. To verify our approach 
for problems of astrophysical interest, we consider the homologous 
expansion phase of SN la ejecta. 

In the expelled material of SNe la, radioactive decays emit 
7-rays, which interact strongly with the ejecta. This coupling af- 
fects the expansion dynamics leading to changes in the density and 
velocity profile that could influence the observable display. The 
question of whether this radiation-hydrodynamical effect is impor- 
tant for theoretical determinat ions of SN light curves and spectra 
has already been addressed bv lPinto & Eastman] (feOQO) using a se- 
ries of analytical estimates. This study concluded, that the radia- 
tive influence on the expansion of the ejecta is marginal and does 
not cause significant chang es in the light curve . These findings 
were verified in the study of IWooslev e t al. (2007), which involved 
detail ed radiation-hydrodyn amical calculations with the STELLA 
code telinnikov et al.ll2006t) . Changes on a 10 per cent level in the 
density stratification of the ejecta were identified when including 
the radiation-matter coupling. Here, we re-address this problem for 
the purpose of testing our radiation-hydrodynamical approach in 
astrophysical environments. 

5.1 Code Extensions 

For the application to supernova ejecta, our program is extended 
to account for the energy release accompanying radioactive de- 
cays. The determination of the energy injection into the radiation 
field involves tracking the abundances of 56 Ni and its daughter nu- 
cleus 56 Co. Both radio-nuclides are proton-rich and decay through 
electron-capture reaction^] down to stable 56 Fe. These decay re- 
actions occur on the characteristic e-folding time scales of inj = 
8.80 d for 56 Ni -> 56 Co and t Co = 113.7 d for 56 Co 56 Fe 
and are accompanied by a cascade of 7-ray emissions. The 7- 
photons interact with the surrounding material through Compton- 
scattering, production of e + e~ pairs and the photoelectric effect. 
As a result of these interaction mechanisms, the 7-radiation heats 

4 Note that the 56 Co dec a y proc eeds via positron emission in about 20 per 
cent of cases. As in|Lucy| |2005j). we assume instantaneous local annihila- 
tion of the positrons yielding a pair of 7-rays. We neglect the kinetic energy 
of the positrons. 



the surrounding material and the energy is re-radiated as quasi- 
thermal emission. Assuming instant thermalisation we simplify the 
7-interaction processes by relying only on one effective a bsorp- 
tive opacity dSufherland & Wheelenll984l ISwartz et al.lll995l) . We 
model the net effect of injecting energy into the thermal radia- 
tion field by a grey Monte Carlo transport step. For each simu- 
lation time step, the number of radioactive decays is determined 
and an adequate number of Monte Carlo packets representi ng the 
emitte d 7-energy is created. For this purpose we follow iLucvl 
(2005) and integrate over the 7-sp ectrum of the decay reactions 
(see lAmbwani & Sutherland]! 19881 table 1), obtaining a total en- 
ergy release in the form of 7-radiation of i?Ni = 1.728 MeV and 
Eco = 3.566 MeV per decaying 56 Ni and 56 Co nucleus respec- 
tively. The 7-packets propagate in the same manner as the Monte 
Carlo packets describing the thermal radiation field, but the inter- 
actions with the ejecta material are described by different opaci- 
ties. Each 7-packet that undergoes an absorption is automatically 
transformed into a thermal radiation packet, which is then treated 
according to the methods laid out in Section [3~!8l 

5.2 Toy Simulations 

To verify the correct operation and the validity of the adopted sim- 
plified treatment of the radioactive energy inje ction, we pre sent the 
results of a simple toy simulation following ILucvl d2005h . In his 
study, Lucy developed a Monte Carlo radiative transfer method to 
determine spectra and light curves in SNe la. In the verification 
process, Lucy performed a grey radiative transfer simulation for a 
simplified ejecta structure under the assumption of homologous ex- 
pansion. A comparison of the bolometric light curve determined in 
the Monte Carlo si mulation with the results of a moments-equation 
solution technique (Castor 1972) yielded excellent agreement. Us- 
ing the results of ILucvl ( 120051) as a reference, we test the imple- 
mentation of the radioactive decay mechanism and the 7-transport 
module in our Monte Carlo radiation hydrodynamics program. In 
addition, this simulation provides yet another test for the accurate 
operation of the Monte Carlo radiative transfer procedures . Partic- 
ularly, the determination of the emergent light curve is sensitive to 
the correct implementation of first-order relativistic effects, such as 
angle aberration and Doppler-shifts, due to the high ejecta veloc- 
ities in this appl ication. The parameters for this test calculations, 
as described in Lucv (2005), are adopted from a model SN pre- 
sented in iPinto & Eastman! ( 2000). We consider a SN with a to- 
tal ejecta mass of M = 1.39M0. The ejecta are in homologous 
expansion with a maximum velocity of it max = 10 4 kms -1 and 
a uniform density of p = 3.79 x 10 -11 gem -3 . In the ejecta, 
the radioactive 56 Ni is assumed to be strongly concentrated in 
the core, resulting in a 56 Ni distribution of fm(M r ) = 1 for 
M r < 0.5 Mq that linearly drops to zero at M r = 0.75 M . 
Here, M r denotes the mass contained within a sphere of radius 
r. Throughout the entire material, a constant interaction opacity 
for the the rmal r adiati on field of x/p = 0.1 cm 2 g -1 is used. 
Following ILucvl d2005l) . we assume radiative equilibrium which 
allows us to simplify the interaction mechanism to only include 
scattering events. The grey absorption cross section for 7-radiation 
of K~,/ p = 0.03 cm 2 g -1 is adopted from [Sutherland & Wheeierl 
dl984l) and lAmbwani & Sutherland] dl988l) . 

We start the Monte Carlo simulation at t — 3.0 d. At much 
earlier times, the high optical thickness of the ejecta material pre- 
vents an efficient use of the Monte Carlo radiative transfer meth- 
ods. We bridge the time between the explosion to the start of the 
simulation by an analytic homologous expansion calculation. Here, 
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Figure 7. Initial conditions for our Monte Carlo simulation of the iLucvl 
model SN ejecta at t = 3.0 d. The time before the onset of the 
simulation was bridged by an analytic homologous expansion calculation, 
resulting in the temperature profile shown as a blue solid line. The dashed 
red line describes the distribution of 56 Ni. 
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Figure 8. Comparison between the bolometric light curve calculated with 
our Monte Carlo program (red) and the results of the Lucy (2005) calcula- 
tion (blue). In addition, the rate at which the 7 radiation deposits its energy 
in the thermal radiation field is illustrated for both calculations (magenta 
and cyan). As a reference, the actual energy release in the decay reactions 
in the form of 7 radiation is also shown as the green line. 



we assume that the entire energy released in the radioactive de- 
cay reactions of 56 Ni and 56 Co immediately thermalises, raising 
the gas temperature to the values shown in Fig. [7] at the time of 
the onset of the Monte Carlo simulation. This figure also illus- 
trates the distribution of radioactive 56 Ni after t he initial hom ol- 
ogous expansion phase. Since the calculations of lLucvl ( 120051) did 
not include any radiation-hydrodynamical coupling, this interac- 
tion mechanism was also switched off in our test simulation. Con- 
sequently, the radiative transfer is not affecting the homologous ex- 
pansion of the ejecta material. 

Fig. [8] presents the bolometric light curve determined in our 
Monte Carlo simul ation and shows the comparison with the re- 
sults of lLucvl J2005h . The excellent agre ement betwe en our simula- 
tion and the calculations performed by iLucvl d2005l) demonstrates 
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Figure 9. The temperature profile (green line) of the W7 model after dis- 
cretization of the original data set onto a spherical grid of 2000 cells and 
performing an homologous expansion calculation up to t = 3d after ex- 
plosion. In addition, the mass fraction of the radioactive 56 Ni is illustrated 
with the dashed blue line. It illustrates the characteristic concentration of 
the 56 Ni in the extended shell in the intermediate ejecta regions of the W7 
model (3 X fO 8 cms -1 < u < 10 9 cms -1 ). All elements with Z > 20 
are denoted as iron group elements (IGE) whose combined distribution is 
shown as the red dashed line. 

the operation of the energy injection process together with the 7- 
transport scheme and once more verifies the accurate performance 
of our Monte Carlo approach as a whole. 

5.3 Application to the W7-Model 

With the SN specific extensions, our program is well suited to ex- 
amine the influence of the radiation-hydrodynamical coupling on 
the observed bolometric light curve. For this purpose we consider a 
S N that is described by the well-known W7 model, first presented 
bv lNomoto et al]! 19841) and extended b y lar ger nuclear networks i n 
the studies of iThielemann et al.l (l986) and llwamoto et al] d 19991) . 
This parameterized one-dimensional SN la explosion model repro- 
duces important obser vables and has thus b ecome a standard ref- 
erence in the literature. tNomoto et al.l (fl984) determined the evolu- 
tion of the fluid dynamical state of the ejecta material until t = 20 s 
after the ignition of the thermonuclear explosion. Characteristic for 
the W7 model is the concentration of the bulk of the 56 Ni in an 
extended shell, spanning the velocity region from 3 x 10 8 cms -1 
to 10 9 cm s _1 , instead of being concentrated in the core. To accu- 
rately resolve this 56 Ni shell, we discretize the original W7 pro- 
file to a one-dimensional spherical grid with 2000 equidistant cells. 
In addition, we slightly adjust the velocity profile of the original 
model to match exactly the homology condition, allowing us to 
make a clear and unperturbed identification of the influence of the 
radiation-matter coup l ing on the ejecta dynamics. As in our cal- 
culation for the lLucvl ( 2005) toy model, we start the Monte Carlo 
radiation hydrodynamics calculation at t = 3d assuming perfect 
homologous expansion behaviour of the ejecta up to this point. The 
resulting temperature profile after this pure homologous phase at 
t — 3d is shown in Fig. [9] together with the distribution of 56 Ni. 
The density stratification at this time is visualised in Fig. [10] In 
analog y to the m odel simulation presented in Section |5T2l we fol- 
low lLucvl 120051) and assume that the ejecta material is in radiative 
equilibrium. Consequently, we can simplify the interaction mecha- 
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Figure 10. The density profile of the W7 model at the start of our Monte 
Carlo simulation at t = 3 d (blue, left scale) and the scattering cross-section 
per gram (red, right scale) obtained from Equation 0511 . 
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Figure 11. Temporal evolution of the mean intensity of the radiation field 
in the radiation-hydrodynamical simulation of the W7 model. Initially, the 
radiation field is confined to the 56 Ni shell but with decreasing optical depth 
it penetrates the inner and outer regions of the ejecta. 



nism between the thermal radiation field and the surrounding ma- 
terial by a pure scattering description. The grey scattering opacity 
is not constant throughout the ejecta radius, but follows the dis- 
tribution of heavy elements, which constitute the dominant opac- 
ity source for the thermal radiation p hot ons. In parti cular, we fol- 
low |Mazzair^Podsiadi^s3 120061) and lSiml l2007t) in setting the 
scattering opacity to 



X/p = iV(0.9XiGE + 0.1), 



(51) 



where the iron group (IGE) involves all elements heavier than cal- 
cium (the combined distribution of these elements can be read off 
from Fig. [9}. The scaling factor N is chosen to ensure a mean in- 
teraction cross section of (x/p) = 0.1cm 2 g _1 . Fig. 1101 shows 
the effect of the heavy elements on the interaction cross section 
throughout the ejecta material. All interactions of the 7-radiation 
are described by a single constant absorption opacity of ft 7 /p = 
0.03 cm 2 g _1 (see Section [5"!2l l. With these parameters, the tempo- 
ral evolution of the hydrodynamical state of the W7 ejecta is calcu- 
lated up to t = 45 d. Fig. QT] illustrates the behaviour of the radia- 
tion field, quantified by the mean intensity J, during this period of 
time. At the beginning, the radiation field is concentrated at the lo- 
cation of the 56 Ni-shell, where initially most of the 7-interactions 
take place due to the high optical depth of the ejecta. As time pro- 
gresses, the ejecta become increasingly transparent to the radiation 
field, which begins to penetrate the inner regions and propagates to 
the outer edge of the supernova ejecta. Here, the Monte Carlo pack- 
ets escape and are recorded to determine the bolometric light curve. 
Due to the radiation pressure, the initially confined radiation field 
leaves imprints on the velocity and density profiles of the ejecta ma- 
terial as it propagates in- and outwards. The changes in the ejecta 
structure are indicated in Fig. [12] which shows the velocity and 
density profiles with respect to a purely homologous expansion that 
would occur in the absence of the radiation-hydrodynamical cou- 
pling. The radiation pressure accelerates the ejecta motion in the 
outer parts of the 56 Ni shell and decelerates the expansion in the 
inner regions. Both effects dilute the 56 Ni shell and pile up material 
at the edges of this region. The influence of the radiation pressure 
is strongest in the early phases of the expansion, where the high 
optical depth causes a strong coupling of the radiation field to the 
surrounding material and stalls as the ejecta become more and more 
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Figure 12. Illustration of the influence of the radiation field on the evolu- 
tion of the fluid state. In the upper panel, the velocity profile in the ejecta 
is displayed with respect to a purely homologous expansion for a variety 
of temporal snapshots (see labels in figure). In the lower panel, the density 
stratification is displayed in a corresponding fashion. The de- and accel- 
erating effects of the radiation pressure are clearly visible in the velocity 
profiles, as is the resulting dilution of the 56 Ni shell in the lower panel. 



transparent. In total, the radiation pressure induces deviations from 
the purely homologous density profile on the order of 10 per cent, 
which are compatibl e with the findings of the previous study by 
IWooslev et all 12007). The structural changes in the density strat- 
ification are, however, not prominent enough to significantly alter 
the shape of the emergent bolometric light curve. Neglecting the 
radiation-hydrodynamical coupling and assuming a purely homol- 
ogous expansion yields nearly the same bolometric light curve as 
a detailed calculation. A comparison of both simulations is shown 
in Fig. [T3] Even if the bolometric light curve remains unaffected 
by changes in the fluid state, colour light curves may be affected 
since the radiation pressure changes the velocity of ejecta regions 
in which different elements are concentrated. However, this effect 
cannot be studied with the current grey implementation of our ap- 
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Figure 13. Comparison between bolometric light curves calculated for the 
W7 model with the radiation-hydrodynamical coupling (blue line) and with 
a pure radiative transfer calculation and homologous expansion of the ejecta 
material (red line). The emergent 7-light curve is also displayed for both 
simulations (green and cyan, respectively). No significant influence of the 
radiation-matter coupling or the resulting changes in the ejecta structure can 
be identified in the bolometric light curves. 



Figure 14. Illustration of the influence of smoothing on the density evolu- 
tion of the ejecta. The density profiles at t = 14.3 d are displayed with re- 
spect to a purely homologous expansion (cf. Fig 1121 . We show the results of 
two simulations: one that employs smoothing over the neighbouring 30 cells 
(red line) and one with very mild smoothing over 3 cells (blue line). As this 
comparison shows, the smoothing approach removes the strong stochastic 
cell-to-cell fluctuations, but retains the overall physical variation. 



proach to radiation hydrodynamics and remains to be re-addressed 
with a chromatic version of the code. 

All simulations presented in this section required about 3.5 x 
10 13 floating point operations. In these calculations the initial radi- 
ation field at t = 3 d was discretized by 10 6 Monte Carlo packets. 
Due to the energy release in the decay reactions and the outflow 
of radiation packets at the outer edge of the ejecta, the number of 
active packets varied greatly during the simulation, but a maximum 
of 2.43 x 10 6 packets were used to represent the radiation field at 
any time. 

Despite the large number of packets, the Monte Carlo estima- 
tors for the radiation force are subject to a significant level of sta- 
tistical noise. In this particular case, the entire energy- momentum 
transfer is determined by the radiative flux in the co-moving frame 
(cf. Eouationsl42land|43l>. As pointed out in Sections [33] and [XTT] 
our discretization approach is not ideally suited to accurately cap- 
ture this quantity in regions where the radiation field is nearly 
isotropic, such as the inner parts of the 56 Ni shell. However, de- 
spite being subject to a considerable level of Monte Carlo noise, 
the radiative flux captures the expected effect of the radiation field 
inflating the 56 Ni bubble (see Figure [T4t. Thus, we employed a 
smoothing kernel (cf. Section [3. Ill in our simulations, averaging 
over 30 neighbouring cells and thereby avoiding difficulties in the 
hydrodynamical solver. In the absence of this smoothing, cell-to- 
cell fluctuations in the fluid state lead to convergence problems in 
the Riemann solver that could not be averted by increasing the num- 
ber of active packets on computationally feasible scales. Using the 
smoothing capability gives a significant reduction of the cell-to- 
cell fluctuations without damaging the physical variations in the 
quantities, as Fig. [14] illustrates. Here, the density stratification at 
t = 14.3 d is shown for both a simulation that includes smoothing 
over 30 neighbouring cells and for one with very mild smoothing 
over 3 cells only, proving that the approach preserves the physical 
signal. 



6 DISCUSSION 

In this paper we presented a Monte Carlo approach to radiation hy- 
drodynamical problems in astrophysical environments. By combin- 
ing Monte Carlo radia tive transfer methods t hat rely on the indivisi- 
ble packet formalism l l Abbott & Lucvl 19851) with the fini te- volume 
hydrodynamical technique PPM dColella & Woodward 1 19841) we 
have aimed to retain the benefits of the Monte Carlo machinery for 
the modelling of complex interaction physics and arbitrary geome- 
tries. Here, our main focus lay on the development and presentation 
of the necessary numerical tools and on demonstrating the oper- 
ation of this method, its physical accuracy and its computa tional 
feasibility. By usin g volume-based Monte Carlo estimators jLucvl 
1 19991 I2002L 120031) in the reconstruction of radiation field charac- 
teristics, the maximum amount of information is extracted from the 
propagation behaviour of Monte Carlo packets and the Monte Carlo 
noise is minimized. 

A series of toy calculations has been performed to test the 
operation of the main components of our Monte Carlo radiation 
hydrodynamical method. In particular, the simulation of radiative 
shocks verified the accuracy of our approach to a standard ra- 
diation hydrodynamical test. As expected, due to the nature of 
the Monte Carlo method, calculations in optically thick environ- 
ments are time-consuming but feasible and accurate, as the radia- 
tive shock examples showed. In general, all calculations were com- 
pleted within hours to a day on a single desktop CPU. However, due 
to the very efficient scaling behaviour of the Monte Carlo algorithm 
to large numbers of processor cores, a future parallel implementa- 
tion of the method provides the scope for significant decreases in 
the run time. 

The application to SN la ejecta successfully demonstrated the 
operation of our method to an astrophysical problem for which this 
method was primarily developed. In this exercise, the influence of 
the radiation-matter coupling on the density stratification during 
the near-homologous expansion phase of the ejecta has been inves- 
tigated. T he results we o b tained are in agreement with the previous 
study by IWooslev et all d2007l) who used the radiation hydrody- 
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namics code STELLA dBlinnikov et alj2 006). The induced changes 
in the ejecta structure, however, were confirmed to have no sig- 
nificant influence on th e bolometric light curve, as predicted by 
IPinto & Eastmanl feOOCh . 

Despite the agreement of our results with previous studies, 
the SN la application also illustrated some difficulties of our ap- 
proach. Our discretization scheme into packets of radiative energy 
allows us to easily construct all relevant radiation field character- 
istics and can be generalized to a fully frequency-dependent trans- 
port treatment in a straight-forward manner. However, in regions 
where the radiation field is close to LTE or to isotropy, the radia- 
tion force components are subject to considerable statistical fluc- 
tuations due to this discretization choice. Here, we suppressed this 
noise component by applying a smoothing kernel. Although be- 
yond the scope of this work, in the future further reduction of the 
statistical noise sh ould be explored by incorp orating implicit Monte 
Carlo techniques jFleck & Cummingslll97ll) or a Monte Carlo ra- 
diative transfer approach that is based on the difference formulation 
dSzoke & Brooksll2005l; iBrooks et aT] |2005). For such schemes, it 
will be important to consider how all the physical processes neces- 
sary to adequately address particular astrophysical applications can 
be implemented. 

As the main aim of this work was to establish the methods 
and the numerical framework of the Monte Carlo radiation hydro- 
dynamics approach, all calculations were performed with a sim- 
plified implementation of the method. In particular, we restricted 
our tests to one-dimensional geometries and the radiative transfer 
was performed in a grey approximation with a very simple opac- 
ity prescription. However, we stress that these simplifications were 
made only to reduce the complexity and the computational effort 
of the test simulations. They do not affect the applicability and 
operation of the Monte Carlo radiation hydrodynamical approach 
itself. In the future, we aim to introduce more sophisticated opac- 
ity prescriptions and make the transition from grey transport to a 
fully frequency-dependent radiative transfer scheme. This will add 
a further level of sophistication to the method, but will not impact 
the operation of the radiation-matter coupling in our approach. The 
tools to realise the frequency-dependent transfer have already been 
developed and are provide d for examp le in the framework of the 
macro-atom formalism by iLucvl d2002l) . In addition to improving 
the physical accuracy of the radiative transfer, our long-term efforts 
will be directed towards a generalised implementation for multi- 
dimensional problems. With this generalization our radiation hy- 
drodynamical method will include all major capabilities that have 
already made the Monte Carlo technique a very successful and re- 
warding approach for pure radiation transport applications. 
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